#ifndef TEST_MATRIX_H
#define TEST_MATRIX_H

#include "TESTS/utilities_for_tests.h"
#include "FEM/CONTAINERS/framedmatrix.h"
#include "FEM/CONTAINERS/simplematrix.h"

#define VIENNACL_WITH_OPENCL
#include "viennacl/matrix.hpp"
#include "viennacl/linalg/prod.hpp"

#include <stdio.h>

#include <QTime>

bool matrixTest1()
{
    /*FramedMatrix<double> m(3,3);
    m[2][2] = 1;
    for(int i=0;i<3;++i)
    {
        for(int j=0;j<3;++j)
        {
            printf("%f ",m[i][j]);
        }
        printf("\n");
    }
    message((char*)"matrixTest1 a",m[2][2] == 1);
    FramedMatrix<int,4> m2(9,9);
    for(int i=0;i<9;++i)
        for(int j=0;j<9;++j)
            m2[i][j]=j+i*9;

    for(int i=0;i<9;++i)
    {
        for(int j=0;j<9;++j)
        {
            printf("%2d ",m2[i][j]);
        }
        printf("\n");
    }
    message((char*)"matrixTest1 b",m2[0][2] == 2 );
    m2.transpose();
    for(int i=0;i<9;++i)
    {
        for(int j=0;j<9;++j)
        {
            printf("%2d ",m2[i][j]);
        }
        printf("\n");
    }
    return message((char*)"matrixTest1 c",m2[2][0] == 2 );*/
    return true;
}

bool matrixTest2()
{
    /*FramedMatrix<double> m1(3,3);
    m1[0][0] = m1[1][1] = m1[2][2] = 2;
    FramedMatrix<double> m2(3,3);
    m2[0][1] = m2[1][2] = m2[2][0] = 1;
    FramedMatrix<double> m3(3,3);
    FramedMatrix<double>::product(m3,m1,m2);
    for(int i=0;i<3;++i)
    {
        for(int j=0;j<3;++j)
        {
            printf("%f ",m3[i][j]);
        }
        printf("\n");
    }
    return message((char*)"matrixTest2 a",m3[1][2] == 2 );*/

    //FramedMatrix<double> m4 = m1 + m2;
    return true;
}

bool matrixTest3()
{
    SimpleMatrix<double> m1(3,3);
    m1[0][0] = m1[1][1] = m1[2][2] = 2;
    SimpleMatrix<double> m2(3,3);
    m2[0][1] = m2[1][2] = m2[2][0] = 1;
    SimpleMatrix<double> m3 = m1 + m2;
    message((char*)"matrixTest3 a",m3[1][2] == 1 );
    SimpleMatrix<double> m4(m1);
    message((char*)"matrixTest3 b",m4[0][0] == 2 );
    SimpleMatrix<double> m5 = m1;
    message((char*)"matrixTest3 c",m5[1][1] == 2 );
    m4 = m1 + m2;
    m5 = m3;
    message((char*)"matrixTest3 d",m4[1][2] == 1 && m5[1][2] == 1);
    m4 = m1 * 2 + 6;
    message((char*)"matrixTest3 e",m4[0][0] == 10 );
    m5 = 2 * m1 * m2;
    message((char*)"matrixTest3 f",m5[1][2] == 4 );
    SimpleMatrix<int> a1(5,2);
    a1[0][0] = a1[1][1] = a1[2][1]= a1[3][0] = a1[4][1]= 2;
    SimpleMatrix<int> b1(2,3);
    b1[0][0] = b1[1][1] = b1[0][2]= 3;
    SimpleMatrix<int> I(2,2);
    I[0][0] = I[1][1] = 1;
    SimpleMatrix<int> c1 = a1 * 4 * I * b1 - 4;
    return message((char*)"matrixTest3 g",
            c1[0][0] == 20 && c1[0][1] == -4 &&
            c1.nCols == 3 && c1.nRows == 5);
}

bool matrixTestSpeedComparison()
{
    QTime t1, t2;
    int N = 512;
    int M = 1024;
    {
        SimpleMatrix<float> smA(N,M);
        SimpleMatrix<float> smB(M,M);
        SimpleMatrix<float> smR(N,M);
        for(int i=0; i<N; i++)
            for(int j=0; j<M; j++)
            {
                smA(i,j)=qrand();
                smB(j,i)=qrand();
            }
        t1.start();
        smR = smA * smB;
        printf("Simple matrix multiplication, time elapsed:     %d ms\n", t1.elapsed());
    }
    {
        viennacl::matrix<float> vmA(N,M);
        viennacl::matrix<float> vmB(M,M);
        viennacl::matrix<float> vmR(N,M);
        for(int i=0; i<N; i++)
            for(int j=0; j<M; j++)
            {
                vmA(i,j)=qrand();
                vmB(j,i)=qrand();
            }
        t2.start();
        vmR = viennacl::linalg::prod(vmA,vmB);
        printf("Viennacl matrix multiplication, time elapsed:   %d ms\n", t2.elapsed());
    }

    return message((char*)"matrixTestSpeedComparison",
            true);
}

#endif // TEST_MATRIX_H
